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Numerical  modeling  of  two-  and  three-dimensional  low  Reynolds  number  gas  flows 
from  small  nozzles  has  been  performed  using  the  direct  simulation  (DSMC)  method. 
The  objective  of  this  effort  is  to  gain  an  improved  understanding  of  performance 
and  plume  interaction  phenomena  for  low  thrust  devices,  and  thus  improve  the 
design  and  optimization  process  for  a  variety  of  micro-propulsion  systems. 
Simulations  were  performed  for  a  wide  range  of  flow  parameters  using  the  SMILE 
parallel  DSMC  code.  Validation  has  been  conducted  through  comparison  of  mass 
flow  and  thrust  values  obtained  numerically  with  results  of  experimental 
measurements  carried  out  recently  by  AFRL  researchers .  Good  agreement  observed 
between  simulation  and  experiment  provided  credibility  to  important  findings 
related  to  plume-plume  and  plume-surface  interactions.  In  particular,  a  large 
effect  of  plume  interaction  on  surface  back  force  and  net  nozzle  thrust  has  been 
shown  for  low  Reynolds  flows.  The  influence  of  the  gas  surface  interaction 
accommodation  coefficient  on  nozzle  performance  was  also  investigated.  An 
example  of  the  plume  flow  from  a  4cm  diameter  nozzle  interacting  with  a  plate  is 
shown  in  Figs.  1  and  2  where  the  pressure  field  and  surface  pressure 
distribution  are  shown  for  a  throat  Reynolds  number  of  60  and  nitrogen 
propellant . 
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Here  Rp  is  the  particle  radius,  ng  is  the  gas  number  density, /(ur)  is  the  gas  velocity  distribution  function  in  a  particle- 
centered  coordinate  frame,  m  is  the  mass  per  gas  molecule,  cr  is  the  relative  speed,  x  is  the  particle  thermal 
accommodation  coefficient,  Tp  is  the  particle  temperature  and  kB  is  Boltzmann’s  constant. 

In  considering  the  Green’s  function  Fp(ur)  for  the  force  imparted  by  the  gas  on  a  rotating  particle,  we  first 
decompose  Fp  into  a  term  which  accounts  for  momentum  transfer  to  the  particle  from  incident  gas  molecules,  a  term 
for  momentum  exchange  during  specular  reflection,  and  a  corresponding  term  for  diffuse  reflection.  It  can  be  shown 
that  the  second  term  is  zero  and  that  only  the  third  term  may  be  influenced  by  particle  rotation.  Next,  using  the 
velocity  superposition  principle,  the  momentum  transfer  due  to  diffuse  reflection  is  separated  into  a  component 
independent  of  rotation  and  a  component  which  depends  on  the  particle  angular  velocity  cop  but  is  independent  of  Tp. 
This  latter  component  is  the  only  contributor  in  a  time-averaged  sense  to  the  force  Green’s  function  Fro,  due  to 
particle  rotation,  and  allows  the  Green’s  function  for  the  total  force  on  the  particle  to  be  computed  as  Fp  =  Fnr  +  Frot. 
Following  the  above  logic,  Frot  is  equal  to  the  product  of  the  collision  frequency  for  diffusely  reflecting  collisions 
and  the  rotational  component  -m<ut>  of  the  average  momentum  transferred  to  the  particle  during  diffuse  reflection. 
Here  ut  is  the  tangential  velocity  of  the  particle  surface  at  the  collision  point  in  an  inertial  particle-centered  reference 
frame.  In  order  to  evaluate  <ut>  we  express  ut  as  a  function  of  the  angle  0  between  -ur  and  the  outward  surface 
normal  unit  vector  n0  at  the  collision  point,  then  multiply  by  the  distribution  function  for  0  and  integrate  over  all 
particle  surface  elements  where  a  collision  may  occur.  We  find: 

9tt 

cq,  xu,  (2) 

By  adding  (1)  and  (2)  and  integrating  over  all  gas  velocity  space,  we  then  determine  the  total  force  on  the  particle. 

In  the  absence  of  any  gas  vorticity  effects  or  forces  not  associated  with  interphase  collisions,  the  rotating  particle 
will  experience  a  damping  moment  about  the  particle  center  which  will  tend  over  time  to  reduce  the  magnitude  of 
(Bp.  Depending  on  the  angle  between  (Dp  and  the  gas  bulk  velocity  relative  to  the  particle,  the  moment  may  also  have 
a  component  orthogonal  to  eop  which  tends  to  change  the  axis  about  which  the  particle  rotates.  This  moment  can  be 
determined  using  the  same  Green’s  function  approach  as  above.  Again,  incident  and  specularly  reflected  gas 
molecules  do  not  contribute  to  the  moment,  and  the  average  moment  contribution  during  diffuse  reflection  depends 
only  on  the  tangential  velocity  of  the  particle  surface.  The  Green’s  function  Mp(ur)  for  this  moment  can  then  be 
expressed  as  the  negative  product  of  the  collision  frequency  for  diffuse  reflection  and  the  average  angular 
momentum  mRp<ncxu,>  about  the  particle  center  imparted  on  diffusely  reflected  gas  molecules.  The  averaging 
operation  is  performed  through  integration  over  the  particle  surface,  to  find  the  following  result: 

H  (>0 = ^l^mn/fu^TC,  ^ u,  -  3  cq,  (3) 

Through  energy  conservation  arguments,  the  Green’s  function  for  the  heat  transfer  rate  to  a  rotating  particle  can 
be  expressed  as 

Q,(Ur)=^n/KM  ^mC?-eki,+fero,-^\oMj  (4) 

where  er0(  is  the  average  rotational  energy  of  incident  gas  molecules,  Aro,  is  the  number  of  rotational  degrees  of 
freedom  for  the  gas,  and  ekin  =  2kBTp+  '/2m<UfUt>  is  the  average  kinetic  energy  of  diffusely  reflected  gas  molecules 
in  an  inertial  particle-fixed  coordinate  frame.  By  integrating  over  the  particle  surface,  we  derive  the  general  formula 

Op(iv)=^n/(ur)xcr  ^ rric^ ~  ~  Aoi jk0Tp  +erot +^mo| Rj|  3-^^- 

where  (q,=j<Bp|.  Note  from  (5)  that  a  greater  value  of  ox,  corresponds  to  an  increase  in  heat  transfer  from  the  gas  to  the 
particle.  This  property  can  be  justified  physically  by  considering  a  single  particle  surface  element  in  a  coordinate 
frame  which  is  fixed  to  the  local  tangential  surface  velocity.  In  this  coordinate  frame  the  average  kinetic  energy  of 
diffusely  reflected  gas  molecules  will  be  independent  of  o^,  while  the  average  kinetic  energy  of  incident  gas 
molecules  tends  to  increase  as  cOp2.  Thus,  the  only  effect  of  particle  rotation  here  is  to  increase  the  interphase  energy 
transfer  associated  with  the  kinetic  energy  of  incident  gas  molecules  which  are  involved  in  diffusely  reflecting 
collisions. 

In  the  special  case  where  the  particle  which  moves  at  the  bulk  velocity  of  an  equilibrium  simple  gas  with 
temperature  Tg,  (5)  may  be  integrated  analytically  over  all  gas  velocity  space  to  provide  a  closed  form  solution  for 
the  net  heat  transfer  rate: 
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q, « = : K  vV27tkB  Vm  (4  +  A^k^-T^mn^  (6) 

If  Tp  >  Tg  as  is  typical  for  high-altitude  rocket  exhaust  plume  flows,  then  the  two  bracketed  terms  in  (6)  are  of 
opposite  sign.  Assuming  the  first  term  is  greater  in  magnitude  than  the  second,  we  find  that  particle  rotation  tends  to 
reduce  the  magnitude  of  the  heat  transfer  rate  and  allows  for  a  more  gradual  reduction  in  Tp  over  time. 


DSMC  IMPLEMENTATION 


The  Green’s  functions  given  above  can  be  used  to  allow  rotating  solid  particles  to  be  modeled  within  a  DSMC 
simulation.  For  simplicity  we  neglect  radiative  heat  transfer  and  assume  the  particle  phase  to  be  dilute  and 
chemically  inert,  so  that  temporal  variation  in  the  properties  of  individual  particles  may  occur  only  through  the 
mechanisms  described  above  or  through  particle-wall  interactions.  Representative  solid  particles  are  tracked  along 
with  DSMC  gas  molecules  through  a  computational  grid,  and  during  each  time  step  the  total  force,  moment,  and 
heat  transfer  rate  for  each  particle  is  determined  by  adding  contributions  from  all  gas  molecules  which  are  assigned 
to  the  same  cell.  The  force,  moment  and  heat  transfer  contributions  to  a  particle  from  a  single  DSMC  gas  molecule 
are  computed  by  evaluating  (1),  (2),  (3)  and  (5)  with  the  quantity  n^/fu,)  replaced  by  the  ratio  of  the  local  gas 
molecule  weighting  factor  to  the  cell  volume.  Once  all  Ng  molecules  in  the  cell  have  been  considered,  the  particle 
velocity  up  and  temperature  Tp  are  altered  as  described  in  Ref.  [8],  and  cbp  is  updated  by  the  product  of  the  net 
moment,  the  local  time  step  and  the  inverse  of  the  particle  moment  of  inertia. 

Due  to  the  large  particle  phase  mass  fractions  typical  of  the  high-altitude  plume  flows  of  interest,  accurate 
simulation  requires  consideration  of  momentum  and  energy  coupling  from  the  particles  to  the  gas.  This  is 
accomplished  here  by  probabilistically  modeling  individual  interphase  collisions  and  modifying  the  velocity  and 
rotational  energy  of  all  gas  molecules  involved.  A  collision  selection  scheme  based  on  the  No  Time  Counter  method 
of  Bird  [12]  is  used  to  determine  which,  if  any,  DSMC  gas  molecules  will  collide  with  each  solid  particle  during  a 
given  time  step.  If  a  particular  gas  molecule  is  found  to  collide  with  the  particle,  then  the  collision  will  involve  either 
diffuse  reflection,  with  probability  t,  or  specular  reflection,  with  probability  1-x.  In  the  latter  case,  the  post-collision 
relative  velocity  ur‘  is  sampled  from  an  isotropic  distribution  such  that  the  relative  speed  cr  of  the  molecule  is 
unchanged  during  the  collision,  and  the  molecule  velocity  is  reassigned  to  equal  up+ur .  In  the  case  of  diffuse 
reflection,  a  more  complicated  procedure  is  required  to  find  ur ,  and  the  gas  molecule  rotational  energy  must  be 
altered  as  well. 

If  a  gas  molecule  is  diffusely  reflected  off  a  nonrotating  solid  particle,  then  the  magnitude  and  direction  of  ur*  are 
statistically  independent,  so  each  can  be  determined  separately.  Following  the  procedure  outlined  in  Ref.  [8],  the 
post-collision  relative  speed  cr  is  sampled  from  the  distribution  function 

/(c;)dc;=2itc;3e.P(-Kc;2)dc;  (7) 


using  the  acceptance-rejection  method  [12],  where  p 


).  The  angle  5  between  the  vectors  -ur  and  ur  is  also 


determined  using  the  acceptance-rejection  method,  by  sampling  from  a  numerical  approximation 

7(8)  —  0.0204256-  0.25 1 55s  +  1.10454  -  1.90353  +  0.493  8  52+  1.2485  (8) 

of  the  distribution  function  for  5s[0,  it].  An  azimuthal  angle  a0  required  to  define  the  direction  of  ur*  in  three 
dimensional  space  is  sampled  from  a  uniform  distribution  over  [0,  27i],  and  the  components  of  ur*  are  calculated 
using  formulas  derived  by  Bird  [12]  for  binary  elastic  collisions: 

<9) 

Cr  CrL  (\+A)  J  Cf|_  «+"[) 

The  absolute  velocity  of  the  gas  molecule  is  then  set  to  the  final  value  up+ur*. 

If  a  diffusely  reflecting  collision  involves  a  rotating  solid  particle,  then  Eq.  8  is  not  valid  and  an  alternate 
procedure  must  be  used.  In  this  case  ur  is  found  by  separately  determining  nonrotational  u„r*  and  rotational  urol* 
components,  then  adding  the  two  together.  To  calculate  unr ,  we  first  compute  the  angle  0  between  -ur  and  the 
outward  normal  unit  vector  nc  at  the  collision  point  on  the  particle  surface  by  setting  sin20  equal  to  a  random  number 
on  [0,  1],  A  corresponding  azimuthal  angle  ai  is  sampled  from  a  uniform  distribution  over  [0,  27t],  and  components 
of  nc  are  then  calculated  as 


Wrtf  +  sirflsu .V  =-  -vrccs9-^«±iW)  ,  1, 

i  Cr|_  ('f+wfT  J  cr 


wcoc9  I  sirffifiyoos  Qj  -urwrsin  <aj)  .  (10) 
(v?+w ]f 
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Next,  we  find  the  angle  C,  between  the  vectors  nc  and  unr*  by  sampling  sin2^  from  a  uniform  distribution  over  [0,  1], 
An  azimuthal  angle  a2  is  sampled  from  a  uniform  distribution  over  [0,  27c],  and  Eq.  7  is  used  to  find  the  magnitude 
of  unr .  The  components  of  unr  are  then  determined  through  the  following  equations: 


v>c;  ,  w;=c;  ^ 


(1-nLr 


^(r^yCCS  <^+i\xr^zsin  ct) 

a-" If 


Finally,  the  rotational  component  of  ur*  is  computed  as  urot*  =  Rp  ropxn0 ,  and  the  gas  molecule  velocity  is  reassigned 
to  the  post-collision  value  up+  unr‘+  urot*. 


PLUME  SIMULATION 

To  evaluate  the  importance  of  particle  rotation  effects,  all  procedures  described  above  are  implemented  in  the 
DSMC  code  MONACO  [10],  and  an  axisymmetric  simulation  is  performed  for  the  nearfield  plume  flow  from  a 
Star-27  solid  rocket  motor  expelling  into  a  vacuum.  The  nozzle  has  an  exit  diameter  of  78  cm  and  a  divergence  half¬ 
angle  of  17.2°.  All  gas  and  particle  properties  along  the  nozzle  exit  plane  are  taken  from  Anfimov  et  al.  [11]  Based 
on  values  provided  by  these  authors,  spherical  A1203  particles  of  diameter  Dp  =  0.15,  0.2,  0.3,  0.5,  1,  2  and  3  pm 
together  constitute  a  total  mass  flow  fraction  of  about  30%,  with  particles  of  each  size  distributed  uniformly  over  the 
nozzle  exit.  The  gas  here  is  a  mixture  of  H2,  N2  and  CO,  with  mole  fractions  of  0.38,  0.31,  and  0.31  respectively. 
Effects  of  radiative  heat  transfer,  heterogeneous  reactions,  gas  chemistry  and  particle  phase  change  are  neglected. 
While  each  of  these  phenomena  may  significantly  alter  flow  properties  of  interest,  we  intend  only  to  assess  the 
importance  of  particle  rotation  effects  on  a  representative  plume  flow,  not  to  calculate  flow  properties  for  a 
particular  case  with  the  greatest  possible  accuracy. 

Following  Vasenin  et  al.  [2],  we  assume  that  coalescing  collisions  between  liquid  droplets  within  the  nozzle  are 
the  dominant  contributors  to  particle  rotation.  The  maximum  allowable  particle  angular  velocity  magnitude  is  then 
determined  by  a  criterion  for  the  centrifugal  breakup  of  liquid  droplets.  A  normalized  angular  velocity 


is  set  to  equal  some  critical  value  Qcri(  for  all  particles  at  the  nozzle  exit,  where  mp  is  the  particle  mass  and  a  is  the 
liquid  surface  tension  of  the  particle  material.  While  the  magnitude  of  cop  will  be  constant  for  all  particles  of  equal 
diameter  at  the  nozzle  exit,  the  initial  direction  of  cop  for  each  particle  is  randomly  selected  on  a  plane  normal  to  the 
axis  of  symmetry.  Salita  [1]  uses  momentum  and  energy  conservation  arguments  to  show  that  the  collision-induced 
rotation  rate  of  spherical  liquid  agglomerates  should  be  limited  by  ncrit=3.31.  Following  these  arguments,  we 
perform  plume  simulations  with  initial  particle  angular  velocities  corresponding  to  0=3.3 1  and  Q=0  in  order  to 
evaluate  potential  upper  bounds  for  the  influence  of  particle  rotation  on  various  flow  properties.  For  both  cases, 
particle-to-gas  coupling  is  considered  through  application  of  the  method  described  above  involving  (10)  and  (1 1),  to 
avoid  any  differences  in  results  caused  by  the  numerical  approximation  in  (8). 

The  simulations  use  a  rectangular  grid  domain,  starting  at  the  nozzle  exit  plane  and  extending  10  m  downstream 
in  the  axial  direction  and  4  m  in  the  radial  direction.  Calculations  are  performed  on  16  1.4  GHz  AMD  Opteron 
processors,  with  each  case  requiring  about  3000  CPU  hours.  As  an  example  of  simulation  results  for  the  D=3.31 
case,  Figure  1  shows  average  angular  velocity  magnitudes  for  particles  of  each  size  as  a  function  of  distance  from 
the  central  axis,  along  a  radial  plane  10  m  downstream  of  the  nozzle  exit.  Figure  2  shows  the  variation  in  average 
particle  temperatures  along  this  same  plane,  for  both  the  Q=3.3 1  and  Q=0  cases.  From  a  comparison  of  TP  values  for 
any  radial  location  and  particle  size,  particle  rotation  is  shown  in  Figure  2  to  have  little  measurable  impact  on 
particle  temperatures  in  the  plume.  However,  for  all  but  the  largest  and  smallest  particles  considered,  rotation  does 
account  for  a  roughly  5  to  10  K  increase  in  average  particle  temperatures  within  about  0.5  m  of  the  central  axis.  This 
trend  is  likely  due  to  the  fact  that  the  gas  density,  and  therefore  the  component  of  the  interphase  energy  transfer  rate 
associated  with  particle  rotation,  tends  to  decrease  with  distance  from  the  axis.  The  trend  is  less  pronounced  for  the 
smallest  particles  because  these  rapidly  lose  the  bulk  of  their  rotational  energy  on  entering  the  grid  domain,  and  for 
the  largest  particles  because  of  their  relatively  low  initial  angular  velocities. 

With  this  one  minor  exception,  we  find  a  general  lack  of  dependence  on  particle  rotation  for  all  flow  properties 
of  interest.  Distributions  of  gas  translational  temperature,  particle  and  gas  number  densities,  maximum  divergence 
angles  for  each  particle  size,  and  various  other  properties  show  differences  between  the  two  cases  which  are 
comparable  to  the  expected  statistical  scatter.  We  can  therefore  conclude  that,  at  least  for  the  plume  simulation 
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considered  here,  particle  rotation  may  be  neglected  with  no  significant  loss  in  overall  accuracy.  The  numerical 
methods  described  above  do  maintain  some  utility  in  making  this  conclusion  possible,  and  in  allowing  for  the 
simulation  of  other  two-phase  rarefied  flows,  such  as  MEMS  flows  involving  frequent  particle-wall  collisions, 
where  rotation  effects  may  still  be  important. 
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FIGURE  1.  Variation  in  particle  angular  velocity  magnitudes  along  a  radial  plane  10  m  downstream  from  the  nozzle  exit. 
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FIGURE  2.  Average  particle  temperature  as  a  function  of  distance  from  the  central  axis. 
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